package ru.usu.gv.utils.math.mds;

import ru.usu.gv.utils.math.Matrix;
import ru.usu.gv.utils.math.Vector;

import java.util.ArrayList;
import java.util.List;

/**
 * @author spupyrev
 * 19.11.2009
 * 
 * Computes the mds by fitting inner products
 */
public class ClassicalScaling
{
	private static final int MAX_ITERATIONS = 50;
	private static final double EPS = 1e-8;

	public Matrix mds(Matrix d, int dimensions)
	{
		return mds(d, dimensions, new ArrayList());
	}

	public Matrix mds(Matrix d, int dimensions, List<Vector> eigenVectors)
	{
		assert (d.getColumnCount() == d.getRowCount());
		for (int i = 0; i < d.getRowCount(); i++)
			for (int j = 0; j < d.getColumnCount(); j++)
			{
				assert (Math.abs(d.get(i, j) - d.get(j, i)) < 1e-6);
				assert (d.get(i, i) == 0);
				assert (d.get(i, j) >= 0);
			}

		Matrix B = new Matrix(d);

		B.squareEntries();

		B.doubleCenter();
		B.scale(-0.5);

		Matrix res = new Matrix(B.getRowCount(), dimensions);
		int uIndex = 0;
		for (int i = 0; i <= dimensions; i++)
		{
			Vector u = new Vector();
			u.initRandomUnitLengthVector(B.getRowCount(), i);
			double lambda = powerIteration(B, u);
			//System.out.println("\u019B_" + (i + 1) + " = " + lambda);

			if (i < dimensions)
			{
				if (lambda >= 0.0)
				{
					eigenVectors.add(u);
					Vector x = new Vector(u);
					x.scale(Math.sqrt(Math.max(lambda, 0)));

					for (int j = 0; j < B.getRowCount(); j++)
						res.set(j, uIndex, x.get(j));
					uIndex++;
				}
				else
					dimensions++;

				for (int j = 0; j < B.getRowCount(); j++)
					for (int k = 0; k < B.getColumnCount(); k++)
						B.set(j, k, B.get(j, k) - lambda * u.get(j) * u.get(k));
			}
		}

		return res;
	}

	private double powerIteration(Matrix B, Vector u)
	{
		double lambda = 0;

		double r = 0;
		double limit = 1.0 - EPS;
		//iterate until convergence but at most 'maxIterations' steps
		for (int i = 0; (i < MAX_ITERATIONS && Math.abs(r) < limit); i++)
		{
			Vector x = B.multiply(u);
			lambda = x.normalize();
			r = u.dot(x);
			u.copy(x);
		}

		return lambda * r;
	}

	public static void main(String argc[])
	{
		Matrix d = new Matrix(4, 4);

		d.set(0, 0, 0);
		d.set(0, 1, 1);
		d.set(0, 2, 2);
		d.set(0, 3, 2);

		d.set(1, 0, 1);
		d.set(1, 1, 0);
		d.set(1, 2, 1);
		d.set(1, 3, 1);

		d.set(2, 0, 2);
		d.set(2, 1, 1);
		d.set(2, 2, 0);
		d.set(2, 3, 1);

		d.set(3, 0, 2);
		d.set(3, 1, 1);
		d.set(3, 2, 1);
		d.set(3, 3, 0);

		Matrix res = new ClassicalScaling().mds(d, 2);
		System.out.println(res);
	}
}
